Radiative coupling of two quantum emitters in arbitrary metallic nanostructures

We propose a general formalism beyond Weisskopf–Wigner approximation to efficiently calculate the coupling matrix element, evolution spectrum and population evolution of two quantum emitters in arbitrary metallic nanostructures. We demonstrate this formalism to investigate the radiative coupling and decay dynamics of two quantum emitters embedded in the two hot spots of three silver nano-spheroids. The vacuum Rabi oscillation in population evolution and the anti-crossing behavior in evolution spectrum show strong radiative coupling is realized in this metallic nanostructure despite its strong plasmon damping. Our formalism can serve as a flexible and efficient calculation tool to investigate the distant coherent interaction in a large variety of metallic nanostructures, and may be further developed to handle the cases for multiple quantum emitters and arbitrary dielectric–metallic hybrid nanostructures.


Scientific Reports
| (2022) 12:6901 | https://doi.org/10.1038/s41598-022-10624-y www.nature.com/scientificreports/ In this paper, we propose a general formalism to efficiently simulate the radiative coupling between two quantum emitters in arbitrary metallic nanostructures. More importantly, based on the radiative coupling, we propose an approach to calculate the evolution spectrum and population evolution of these two quantum emitters. As an illustrating application, we investigate the radiative coupling and decay dynamics of two quantum emitters embedded individually in the two nanogaps (hot spots) of three silver nano-spheroids. Due to the enormously enhanced electric field in the two hot spots, despite the strong plasmon damping, the strong radiative coupling between the two quantum emitters can still be manifested by the vacuum Rabi oscillation in population evolution and the anti-crossing behavior in evolution spectrum. This formalism may serve as a flexible and efficient theoretical tool for the distant coherent interaction between two quantum emitters in a large variety of metallic nanostructures, and may be further developed to handle the cases for multiple quantum emitters and arbitrary dielectric-metallic hybrid nanostructures.

Results
Theory. In this section, we deduce a formalism to calculate the temporal evolution of the upper-level-probability amplitudes of two two-level quantum emitters (denoted as A and B ) inside arbitrary metallic nanostructure. By adopting the dipole approximation and the rotating-wave approximation, the Hamiltonian of this system can be expressed as 37,42,43 Here, the first term represents the field energy of the environment in the presence of the metallic nanostructure.f (r, ω) and f † (r, ω) are the bosonic fields that represent the elementary (energy) excitations of the electromagnetic field both in the environment and the metallic nanostructure 44 .
The second term in the Hamiltonian represents the energy of the two quantum emitters. ω i and d i = d idi are the transition frequency and transition dipole moment with magnitude d i and direction d i , respectively, between the excited state |e i � and the ground state g i of the ith ( i = A, B ) quantum emitter located at r i .
The third term in the Hamiltonian represents the interaction between the quantum emitters and the field excitations. σ † i and σ i are Pauli operators of the ith quantum emitter. The electric field operator is separated into two parts as Ê (r) is the electric field operator in the frequency domain 44 . The positive frequency part Ê (+) (r) can be expressed as 42 Here, ε I (r, ω) is the imaginary part of the complex relative permittivity ε(r, ω) of the metallic nanostructure.G(r, r ′ , ω) is the classical Green function (tensor), describing the system response at r to a point source at r ′ , and satisfying the equation together with the boundary condition at infinity. δ(r) is the dyadic δ function. We denote the system states as |a� = e A , g B , 0 , |b� = g A , e B , 0 and |m(r, ω)� = g A , g B , 1 m (r, ω) , with only one excitation at quantum emitter A , B or the bosonic field f m (r, ω) ( m denote the x , y or z component), respectively. The system state evolves as Here, C a (t) , C b (t) and C m (r, ω, t) are the probability amplitudes of |a� , |b� and |m(r, ω)� , respectively. By substituting Eqs. (1), (2) and (4) into Schrödinger equation i ∂ ∂t |�(t)� =Ĥ|�(t)� , and considering the orthonormal relationship among |a� , |b� and |m(r, ω)� , we can obtain (2) C m (r, ω, t) = −iωC m (r, ω, t)+ 1 πε 0 www.nature.com/scientificreports/ We assume that initially there is only one excitation at the two quantum emitters and no excitation at the bosonic fields, i.e., the probability amplitudes for the initial state of the system are We take the Laplace transform (forward and backward Fourier transform) to transform Eqs. (5)-(7) into algebraic equations. The forward and backward Fourier transform of any time-dependent variable C(t) , e.g., C a (t) , C b (t) and C m (r, ω, t) , can be defined, respectively, as

Scientific Reports
Here, the complex frequencies � ± = � ± iη are assumed to contain a real frequency and an infinitely small positive and negative imaginary part ±iη(η → 0 + ), respectively, so that the transform is well-defined 45 .
As derived in the "Methods" section, the forward and backward Fourier transform of C a (t) and C b (t) can be calculated as Here, as derived in the "Methods" section and discussed in detail later, the coupling matrix element is Here, G / * denote G and G * , corresponding to W + ij and W − ij , respectively. Furthermore, the evolution spectrum of the two quantum emitters can be obtained as Finally, the temporal evolution of the two quantum emitters can be calculated via their evolution spectrum as Obviously, the key to investigate the decay dynamics of two quantum emitters in arbitrary metallic nanostructure is the calculation of the coupling matrix element W ± ij (�) , which can be obtained via Green tensor G(r i , r j , �) according to Eq. (13). G(r i , r j , �) describes the system response at r i (the ith quantum emitter) to a point source at r j (the jth quantum emitter) and characterizes the electromagnetic environment of the two quantum emitters, which determines the radiative coupling between them.
For specific metallic nanostructures with simple geometries, e.g., spherical and spheroidal geometries 46 , G(r i , r j , �) can be obtained analytically. However, for arbitrary metallic nanostructures with arbitrary geometries and components, numerical calculations should be adopted to obtain G(r i , r j , �) . Actually, the Green tensor G(r i , r j , �) and hence the coupling matrix element W ± ij (�) can be obtained by simulating the electric fields induced by two point-dipoles, individually.
From Maxwell equations, the electric field induced by an oscillating point-dipole (related to the jth quantum emitter) located at r j with unit magnitude, direction d j and frequency is 47,48 Specially, the electric field component at location r i and along direction d i (related to the ith quantum emitter) is Here, E / * j denote E j and E * j , corresponding to W + ij and W − ij , respectively. This relationship provides a flexible and efficient approach to calculate the coupling matrix element W ± ij (�) , based on the numerical simulation of the electric fields E j (r i , �) induced by two individual point-dipoles.
E j (r i , �) in arbitrary metallic nanostructure can be simulated directly in frequency domain, e.g., by the COM-SOL Multiphysics. Alternatively, at first, E j (r i , t) in time domain can be simulated, e.g., via the finite-difference time-domain (FDTD) method 49 , and then E j (r i , �) can be obtained by Fourier transform or Padé approximation with Baker's algorithm 50 .
2 characterizes the radiative coupling between the ith and the jth quantum emitters, mediated by the electromagnetic field (mainly the plasmon) in the metallic nanostructures. The real and imaginary part of W ± ij (�) corresponds to the collective level shift � ij (�) and the transfer rate Ŵ ij (�) , respectively 26,32,51 Besides, according to Eqs. (13) and (20), we can separately calculate � ij (�) and Ŵ ij (�) via the real and imaginary part of E j (r i , �) , respectively, as 2 characterizes the local coupling between the ith quantum emitter and the electromagnetic field (mainly the plasmon) in the metallic nanostructure. The real and imaginary part of W ± ii (�) corresponds to the level shift � ii (�) and the local coupling strength Ŵ ii (�) of the ith quantum emitter, respectively. Unfortunately, W ± ii (�) can't be directly calculated via Eq. (20). Instead, at first, its imaginary part can be calculated via , then its real part can be calculated via the , via Eqs. (11) and (12), we can successively obtain the evolution spectrum of (14) and (15) and the temporal evolution of Eqs. (16) and (17).
In plasmonic systems, the dominant mechanism is surface plasmon exchange, i.e., excitation of a virtual surface plasmon in the metallic nanostructure by an excited quantum emitter followed by its absorption by the other quantum emitter, rather than direct radiative coupling 31 . The simulated electric field E j (r i , �) corresponds to the excitation and absorption of a virtual surface plasmon. Actually, all system responses including the dominating surface plasmon and other minor electromagnetic responses are totally incorporated in the simulated electric field E j (r i , �) . So this formalism is accurate in calculating the radiative coupling between two quantum emitters.
This formalism is flexible and efficient. It can simulate the radiative coupling and decay dynamics of two quantum emitters at any location, along any polarization direction, with any transition frequency, inside arbitrary metallic nanostructure including, but not limited to, plasmonic cavity and plasmonic waveguide. Besides, only the electric fields at the locations of the two quantum emitters need to be simulated, rather than those of the whole space. This saves massive computational time and memory.

Simulation.
To demonstrate and verify this formalism, we apply it to investigate the radiative coupling and decay dynamics of two quantum emitters in the metallic nanostructure composed of three silver nano-spheroids optimized in Ref. 20 . As shown in Fig. 1, three identical silver nano-spheroids are lined along their elongated axis (along x direction) with axis length a = 13.3 nm . The other two orthogonal axes (along y and z directions) have the same axis length b = 8 nm . The gap widths between any two-neighboring nano-spheroids along the x direction is w = 2 nm . The center of the middle nano-spheroid locates at the origin of the coordinate system. The three nano-spheroids are embedded inside lossless host medium of dielectric constant ε h = 2.2 . Here the refractive index parameters of Ag are obtained by interpolation from the original experimental data in the literature 54 . At the resonant wavelength r = 525.6 nm (resonant frequency r = 2.3589 eV ) of this metallic nanostructure, the electric field are strongly confined and enhanced inside the two gaps, forming two hot spots of electric field, as shown in Fig. 1.
To achieve strong radiative coupling, according to Eq. (20), the two quantum emitters should be individually positioned at the hot spots of electric field. In this three nano-spheroid structure, a natural choice for the locations of the two quantum emitters is the two gap centers, w/2 away from the neighboring nano-spheroids. we take the transition dipole moment of the two quantum emitters as d A = d B = 4.167 × 10 −29 C m , corresponding to the lifetimes in the homogeneous host medium is τ A = τ B = 2 ns . The transition dipole moment of the two quantum emitters are both polarized along x direction to obtain the strongest radiative coupling.
Scientific Reports | (2022) 12:6901 | https://doi.org/10.1038/s41598-022-10624-y www.nature.com/scientificreports/ Based on the electric fields induced by two individual dipole (with unit magnitude, at the locations and along the polarizations of the two quantum emitters, with different frequency), via Eq. (20) and principal value integral, we calculate the coupling matrix element W ± ij (�) between quantum emitter A and B , as shown in Fig. 2. For any i and j , W − ij (�) is the complex conjugate of W + ij (�) . Furthermore, since the three nano-spheroid structure is symmetrical and the two quantum emitters have the same transition dipole moments, we can obtain W ± AA (�) = W ± BB (�) and W ± AB (�) = W ± BA (�). For the imaginary parts related to the local coupling strength Ŵ AA (�) in Fig. 2a and to the transfer rate Ŵ AB (�) in Fig. 2b, there are both two resonant peaks (dips) at 525.6 nm and 477.45 nm . For resonant wavelength of 525.6 nm and 477.45 nm , the electric field distribution is symmetric and antisymmetric about the x = 0 plane, respectively. We focus on the symmetric mode of resonant wavelength r = 525.6 nm , referred to as cavity mode, as shown in Fig. 1. The linewidth for its resonant peak is 10.22 nm (cavity leakage κ = 45.9 meV).
In this paper, we focus on the ideal case that the two quantum emitters have the same transition wavelength A = B . Besides, we assume that at initial time t = 0 , quantum emitter A is in excited state, quantum emitter B is in ground state, and there is no excitation at the bosonic fields, i.e., the initial state of the system is C a (0) = 1 , C b (0) = 0 , C m (r, ω, 0) = 0.
At first, we consider the resonance case that the two quantum emitters are both resonant with the three nano-spheroid structure, i.e., A = B = r = 525.6 nm . By adopting Eqs. (11), (12), (14) and (15), we calculate the evolution spectrum of the two quantum emitters, as shown in Fig. 3a. Since C a (0) and C b (0) are both real, we conclude from Eqs. (11) and (12) that c + a (� + ) and c + b (� + ) are complex conjugate of c − a (� − ) and c − b (� − ) , respectively. In this case, according to Eqs. (14) and (15), the evolution spectrum of the two quantum emitters c a (�) and c b (�) are both real.  www.nature.com/scientificreports/ As shown in Fig. 3a, there are three peaks (dips) in the evolution spectrum c a (�) and c b (�) , respectively. The central peak (dip) at wavelength 528.5 nm (near A = B = 525.6 nm ) is the dark mode, i.e., the antisymmetric state of the two quantum emitters. The central peak is relatively narrow since the linewidths of the two quantum emitters are negligible. In contrast, the two side peaks (corresponding to the upper and lower polaritons) at 508.95 nm and 543.35 nm , with splitting of 34.4 nm ( 154 meV ), are the bright modes composed of the symmetric state of the two quantum emitters and the cavity mode. The linewidths of the two side peaks are relatively large due to the large linewidth of the cavity mode. The large splitting between the two side peaks results from the strong radiative coupling between the two quantum emitters, mediated by the plasmon of this three nanospheroids structure with extremely strong electric field in the hot spots.
The strong radiative coupling is most apparent in time domain. Based on the evolution spectrum c a (�) and c b (�) , via Eqs. (16) and (17), we further obtain the time-dependent probabilities of excitation, i.e., the populations of quantum emitter A and B as P a (t) = |C a (t)| 2 and P b (t) = |C b (t)| 2 . As shown in Fig. 3b, initially, quantum emitter A is in excited state and quantum emitter B is in ground state. For t > 0 , due to the strong radiative coupling between the two quantum emitters mediated by the plasmon (cavity mode), the populations of the two quantum emitters both oscillate at the frequency of g = 80.2 meV , about half of the splitting ( 154 meV ) between the two side peaks in Fig. 3a. The excitation energy is coherently transferred between the two quantum emitters. Meanwhile, due to the cavity leakage κ = 45.9 meV , the amplitudes of the two oscillating populations damp quickly to zero in picosecond scale. The excitation energy finally dissipates mainly due to the ohmic loss of the metallic nanostructure. This pronounced population oscillation is a clear signature of strong radiative coupling, which occurs in a regime where the radiative coupling g is stronger than the plasmon damping κ.  www.nature.com/scientificreports/ Now we turn to investigate the off-resonance case by simultaneously tune the transition wavelength of the two quantum emitters (keeping A = B and r = 525.6 nm ), and calculate the evolution spectrum c a (�) and c b (�) for different emitter-cavity detuning A − r .
As shown in Fig. 4, for positively large emitter-cavity detuning A − r , there are three peaks in both c a (�) and c b (�) . The central peak at A is dark mode, i.e., the antisymmetric state of the two quantum emitters. The side peak approaching A can be mainly attributed to the quantum emitter. The other side peak approaching r can be mainly attributed to the cavity mode. The attribution of the three peaks can also be verified by their linewidths comparing with those of the quantum emitters and cavity mode.
As emitter-cavity detuning A − r decreases, the two side peaks gradually repel each other and can both be attributed to the two quantum emitters and the cavity mode. This behavior is quite similar to the strong coupling system composed of a single quantum emitter and metallic nanostructure 55 .
For zero emitter-cavity detuning, i.e., A = B = r = 525.6 nm , the two side peaks forms two polaritonic states, which is the case of Fig. 3a.
For negatively large emitter-cavity detuning A − r , the evolution spectrum c a (�) and c b (�) are both interfered by the other resonant peak at 477.45 nm , i.e., the antisymmetric mode in Fig. 2. It leads to an extra minor anti-crossing behavior.

Discussion
In this paper, we have proposed a general formalism to calculate the plasmon-mediated radiative coupling between two quantum emitters in arbitrary metallic nanostructure. The coupling matrix element W ± ij (�) can be flexibly and efficiently calculated by simulating the electric fields E j (r i , �) induced by two point-dipoles, individually. Based on the coupling matrix element W ± ij (�) , the evolution spectrum and population evolution of the two quantum emitters can be obtained.
We have demonstrated this formalism to investigate the radiative coupling and decay dynamics of two quantum emitters located in the two hot spots of three silver nano-spheroids. The vacuum Rabi oscillation in population evolution and the anti-crossing behavior in evolution spectrum are clearly observed for both quantum emitters. Obviously, despite the strong plasmon damping, the strong radiative coupling between the two quantum emitters can still be realized in this metallic nanostructure, due to the enormously enhanced electric field in the two hot spots.
This formalism can serve as a flexible and efficient calculation tool to investigate the distant coherent interaction between two quantum emitters in a large variety of metallic nanostructures including, but not limited to, plasmonic cavities and plasmonic waveguides. Besides, it can be further developed to simulate the cases for multiple quantum emitters, which is essential for multiqubit manipulation. It can also be further developed to simulate the case for other nanostructures with arbitrary geometries and components, e.g., the hybrid nanostructures composed of both dielectric and metallic components 56 , where the strong radiative coupling and long-time coherence might be simultaneously realized for quantum information processing. (t)e i� + t dt = −C(0) − i� + c + (� + ).

Methods
. (�+iη−ω)πε 0 c 2 d i · Im[G(r i , r j , ω)] · d j as coupling matrix element between state |i� and state j . If i = j , it means local coupling. If i = j , it means radiative coupling. The detailed derivation of the results can be found from Eq. (30).
Coupling matrix element. The coupling matrix element can be derived as Here, P denotes the principal value integral, and we adopt the following definitions 26,51 � ij (�) can also be expressed as 32 Regarding Eq. (33), Eq. (30) can be further expressed as Here, G / * denote G and G * , corresponding to W + ij and W − ij , respectively. Combining Eqs. (30) and (34), we can derive Eq. (13).

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Code availability
The codes that support the findings of this study are available from the corresponding author upon reasonable request.